Linear Models

Matthew Talluto

03.05.2022

Regression problems

Many applied problems take the form of a set of observations $y$ along with one or more covariates $x$
We want to know how $y$ changes (on average) with $x$, and how much variance there is in $y$

Regression problems

We can draw a line through the points that minimizes the sum of squared residuals $\sum_i (y_i - \mathbb{E}[y|x_i])^2$
$\mathbb{E}[y|x_i]$ is the expectation of $y$ conditional on x
- what do we expect $y$ to be, on average, for any value of $x$

Regression problems

We can draw a line through the points that minimizes the sum of squared residuals $\sum_i (y_i - \mathbb{E}[y|x_i])^2$
$\mathbb{E}[y|x_i]$ is the expectation of $y$ conditional on x
- what do we expect $y$ to be, on average, for any value of $x$
This line has two parameters, an intercept $a$ and a slope $b$
- $\mathbb{E}(y|x_i) = a + bx_i$

Regression problems

We can draw a line through the points that minimizes the sum of squared residuals $\sum_i (y_i - \mathbb{E}[y|x_i])^2$
$\mathbb{E}[y|x_i]$ is the expectation of $y$ conditional on x
- what do we expect $y$ to be, on average, for any value of $x$
This line has two parameters, an intercept $a$ and a slope $b$
- $\mathbb{E}(y|x_i) = a + bx_i$
$y \sim \mathcal{N}(a + bx, s)$
- “y is distributed as”

Regression problems

We can draw a line through the points that minimizes the sum of squared residuals $\sum_i (y_i - \mathbb{E}[y|x_i])^2$
$\mathbb{E}[y|x_i]$ is the expectation of $y$ conditional on x
- what do we expect $y$ to be, on average, for any value of $x$
This line has two parameters, an intercept $a$ and a slope $b$
- $\mathbb{E}(y|x_i) = a + bx_i$
$y \sim \mathcal{N}(a + bx, s)$
- “y is distributed as”
Key assumptions:
- $y$ is iid
- $y$ has constant variance with respect to $x$
- the residuals are normally distributed (or $y|x$ is normally distributed)

Graphing the model

It is always a good idea to graph your model
The graph represents knowns and unknowns as nodes, with relationships as edges
Variables (nodes) may be either random or deterministic
Your model must estimate all unknowns
- random unknowns are usually parameters
- deterministic unknowns must be computed somewhere in your model
- random knowns usually appear as the first parameter is a probability statement like dnorm

Regression likelihood

Your model must estimate all unknowns
- random unknowns are usually parameters
- deterministic unknowns must be computed somewhere in your model
- random knowns usually appear as the first parameter is a probability statement like dnorm

#' params: named vector of parameters
#' data: list or data frame of all data
log_liklihood = function(params, data) {
  # unpack params and data
  a = params['a']
  b = params['b']
  s = params['s']
  x = data[['x']]
  y = data[['y']]
  
  # compute mu, the expectation of y
  mu = a + b*x
  
  # liklihood of y|x,a,b,s
  sum(dnorm(y, mu, s, log=TRUE))
}

Regression prior and posterior

All random variables must appear in a probability statement, i.e. d****()
- knowns belong in the likelihood
- unknowns need a prior

#' params: named vector of parameters
log_prior = function(params) {
  ## here the prior hyperparameters are hard-coded
  dnorm(params['a'], 0, 10, log=TRUE) + 
    dnorm(params['b'], 0, 5, log=TRUE) + 
    dexp(params['s'], 0.2, log=TRUE)  
}

#' params: named vector of parameters
#' data: list or data frame of all data
log_posterior = function(params, data) {
   # s must be positive; if we try an invalid value
   # we have to return something sensible
   # probability is 0, so log probability is -Inf
  if(params['s'] <= 0)
    return(-Inf) 

  log_liklihood(params, data) + log_prior(params)
}

Estimation

data('iris')
iris = iris[iris$Species != "setosa",]
data = with(iris, data.frame(x = Sepal.Length, y = Petal.Length))

param_init = c(a=0, b=0, s=1)
mod_map = optim(param_init, log_posterior, data = data, 
                method="Nelder-Mead", control=list(fnscale=-1))
mod_lm = lm(y~x, data=data)

mod_map$par
##       a       b       s 
## -1.5516  1.0312  0.4598
coef(mod_lm)
## (Intercept)           x 
##      -1.556       1.032
sd(mod_lm$residuals)
## [1] 0.4623

Estimating the entire posterior

Reminder: we want to estimate the joint posterior $pr(a,b,s|x,y)$
Often, the conditional posterior—e.g., $pr(a|b,s,x,y)$—is approximatly normal, and the log-posterior is a parabola
Constrained parameters can often be easily transformed: $pr(\log(s)|a,b,x,y)$ will be approximately normal
It follows that the joint posterior will often follow a multivariate normal
This is a consequence of the central limit theorem

Laplace Approximation

With LA, we find the maximum of the (joint) posterior distribution, and then estimate the Hessian matrix
Hessian is a matrix of second derivatives, tells you how curvy the function is in each direction
The inverse of the Hessian gives us the variance-covariance matrix of a multivariate normal distribution
The posterior is now fully characterized
- fit$par gives us the mean/best-fit for each parameter
- vcv tells us the variance-covariance matrix of the multivariate normal distribution
Now what??

log_posterior_la = function(params, data) {
  ## s must be positive, but optim doesn't know this.
  ## We also want to optimize a parameter that is normal 
  ## working with log(s) can solve both problems
  params['s'] = exp(params['log_s'])

  log_liklihood(params, data) + log_prior(params)
}
param_init = c(a=0, b=0, log_s=0)
fit = optim(param_init, log_posterior_la, data = data, 
            method="Nelder-Mead", 
            control=list(fnscale=-1), hessian = TRUE)
vcv = solve(-fit$hessian)

fit$par
##       a       b   log_s 
## -1.5508  1.0311 -0.7769
vcv
##                a          b      log_s
## a      1.923e-01 -3.038e-02  4.855e-05
## b     -3.038e-02  4.851e-03 -7.420e-06
## log_s  4.855e-05 -7.420e-06  4.994e-03

sds = sqrt(diag(vcv))
sds
##       a       b   log_s 
## 0.43856 0.06965 0.07067